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STATISTICAL  MODELING  OF  THE  AIR-BLAST  ATOMIZATION  IN  THE 
LAGRANGIAN  COMPUTATION  OF  LIQUID  SPRAY 


M.GOROKHO  V SKI 

CORIA  UMR/6614  CNRS  University  of  Rouen,  Mont  Saint  Aignan,  France 
I.  INTRODUCTION 

When  injection  of  liquid  jet  takes  place  into  coflowing  motion  of  high  relative  velocity  gas,  a 
wide  range  of  turbulent  eddies  may  impact  on  the  liquid  jet  causing  its  breakup.  This  regime 
of  atomization  is  often  referred  to  as  the  air-blast  atomization  and  is  widely  used  in  practical 
systems.  The  physics  of  air-blast  atomization  is  very  complex  . In  addition  to  the  gas 
turbulence-induced  breakup,  many  other’s  random  processes  such  as  multiple  droplets 
collision,  turbulence  in  liquid,  variations  in  the  cavitating  flow  inside  the  injector,  etc., 
contribute  to  the  phenomenon  of  breakup.  This  implies  that  at  each  spray  location,  the  spectra 
of  size  of  produced  droplets  can  be  very  large.  Then  the  question  of  what  are  the  probabilities 
of  sizes  that  are  involved  into  atomization,  arises  in  the  breakup  modeling.  Due  to  the 
complexity  of  phenomenon,  it  is  too  difficult  to  disclose  clearly  a dominant  mechanism  of 
air-blast  breakup  with  expectation  of  a characteristic  size  of  droplet.  To  this  end,  the  basic 
idea  behind  the  simulation  undertaken  in  this  paper,  is  as  follows.  The  process  of  air-blast 
breakup  is  considered  in  the  framework  of  cascade  of  uncorrelated  breakage  events  in  series  , 
independently  from  the  initial  distribution  of  sizes.  The  stochastic  modeling  of  droplets 
production  under  this  hypotheses  down  to  the  critical  (or  maximum  stable)  size  is  the  subject 
of  this  paper. 

The  cascade  idea  of  breakup  comes  from  the  early  work  of  Kolmogorov  written  in  1941  6.  In 
this  work,  Kolmogorov  described  the  breakup  of  solid  particles  as  a discrete  random  process, 
where  the  probability  to  break  each  parent  particle  on  a given  number  of  parts  is  independent 
of  the  parent  particle  size.  From  Lyapunov’s  theorem,  Kolmogorov  has  pointed  out  that  such 
a general  assumption  leads  to  the  log-normal  distribution  of  particle  size  in  the  long-time 
limit.  In  this  paper,  the  Kolmogorov’s  discrete  model  has  been  reproduced  in  the  form  of 
evolution  equation  for  distribution  function.  The  asymptotic  solution  of  this  equation  has  been 
applied  to  simulate  the  drop  breakup  alongside  with  Lagrangian  model  of  spray  dynamics. 
Performed  computations  of  air-blast  atomization  are  related  to  a spray  close  to  the  rocket 
engine  configuration. 


II.  KOLMOGOROV’S  (1941)  THEORY  OF  THE  PARTICLE  BREAKUP. 


Let  us  consider  an  ensemble  of  breaking  solid  particles  at  discrete  time  moments  r =0, 1, 2, 

These  time  moments  are  scaled  by  the  breakup  frequency  v (t=vt).  According  to 
Komogorov6,  the  number  of  particles  N(r,t)  of  size  p<r  was  selected  amongst  all  particles 
N(  t ) at  a given  moment  t . The  expectations  of  total  number  of  particles  and  of  particles  of 
size p<  r were  denoted  as  N(t)  and  N(r,t ) correspondingly.  Considering  an  outcome  of 
breakup  per  unit  time  [/,  / + 1]  of  a given  parent  particle  of  size  r , the  mean  number  Q(a ) of 
secondary  particles  of  size  p<ar  ( 0 < a < 1 ) was  introduced.  According  to  hypotheses  of 
Kolmogorov,  the  probability  to  break  each  parent  particle  on  a given  number  of  parts  is 
independent  of  the  parent  particle  size.  In  other  words,  Q{a ) does  not  depend  of  prehistory  of 
breakup  and  is  not  influenced  by  others  parent  particles.  By  this  assumption,  Kolmogorov 
writes: 
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(1) 


Introducing  the  logarithm  of  particle-size  x = In  r , Kolmogorov  pointed  out  that 


T(x,t)= 


jvfev)  ffklA 

N(l)  N{t ) 


(2) 


Further,  denoting  £ = ln«  and  Q(a)  = Q(l)- S(£),  equation  (1)  is  rewritten  by  Kolmogorov  in 
the  following  form: 
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T(x,t  + l)=  \T{x-$,t)dS(4) 


(3) 


By  Lyapunov’s  theorem,  Kolmogorov  stated  that  from  discrete  model  (3),  the  long-time  limit 
form  of  T(.v,/)  tends  to  Gaussian  function.  Then  the  main  result  of  Kolmogorov’s  work  is  that 
N(r,t))  is  asymptotically  governed  by  log-normal  law. 


III.  THE  ASYMPTOTIC  DIFFERENTIAL  FORM  OF  THE  DISCRETE 
KOLMOGOROV’S  MODEL. 


Here  the  discrete  model  (3)  is  represented  by  its  differential  approximation  in  the  long  time 
limit.  Using  parabolic  scaling  of  variables  z = £2 1 , y = £x,  where  £ is  a scaling  parameter 
and  t is  scaled  by  breakup  frequency,  the  equation  (3)  can  be  written  as 


r(y,r+f2)=  \T{y-£Z,r)s(Z)d% 


(4) 


Expanding  both  the  left-hand  side  and  the  expression  under  integral  in  (4),  one  gets 


T(y,T  + £2)=T{y,r)+£2 
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Substituting  these  expansions  in  (4)  and  coming  back  to  variables  tandx , one  yields: 


jgfed , e(c4)-  i 1 

5 1 5 x 2!  3 a 
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where  (£)  = j^s(^)d^  and  |<^2  s(£)d£  are  two  first  moments  of  g . Assuming  that 


i 

the  integral  jjln  a j3  dQ(a)  is  limited,  the  equation  (5)  can  be  written  in  the  long-time  limit 
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£ — > 0 (t  —>o o ),  as 
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dt 
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The  dimensional  time  has  been  used  in  (6).  The  solution  of  (6)  is  Gaussian  function.  This 
repeats  the  main  result  of  Kolmogorov6.  At  the  same  time,  an  influence  of  the  initial 
distribution  before  breakup  starts  can  be  taken  into  account  by  using  (6).  The  solution  of  (6) 
verifies  to  be: 


T’ofa  -{£)vt)dx{) 


(7) 


where  ro(.v„)  is  the  initial  distribution  of  the  logarithm  of  droplet  radius  and  .t0  is  logarithm  of 
radius  of  the  parent  drop.  One  can  rewrite  equation  (6)  for  the  normalized  distribution  of 
radius  /(r): 


(S) 

The  solution  of  this  equation  has  the  following  form: 


(9) 


where  /0(r0)  is  the  initial  distribution  of  droplet  radius  before  breakup  starts. 

IV.  IMPLEMENTATION  INTO  COMPUTATIONAL  CODE  KIVA  II 
IV.  1 General  procedure 

The  right  hand  side  of  equation  (8)  can  be  added  to  the  transport  spray  equation7  as  a source 
term  due  to  drop  breakup.  The  modeling  of  the  spray  equation  is  often  based  on  Lagrangian 
formulation8.  The  spray  is  considered  to  be  composed  of  discrete  parcels  of  particles,  each  of 
which  represents  a group  of  droplets  of  similar  size,  velocity  and  position.  These  groups  of 
droplets  are  followed  as  they  interact  and  exchange  momentum  and  energy  with  surrounding 
gas.  The  basic  ideas  of  this  method,  including  the  modeling  of  turbulent  dispersion  of 


particles,  are  presented  in  9.  Here,  the  Lagrangian  tracking  is  coupled  with  stochastic 
computing  of  breakup.  Two  additional  physical  processes  were  included  in  the  Monte  Carlo 
procedure.  Namely,  the  product  droplet  velocity  has  been  modeled  and  the  breakup  has  been 
considered  down  to  the  local  magnitude  of  the  critical  (or  maximum  stable)  radius,  rcr . The 

liquid  fuel  was  injected  in  the  axial  nozzle  direction  in  form  of  discrete  parcels  of  drops  with 
characteristic  size  equal  to  the  exit  nozzle  radius.  The  injection  velocity  was  determined  from 
the  known  liquid  injection  rate. 

Let  us  consider  the  motion  of  a given  j — th  primary  parcel  that  undergoes  breakup.  Before 
breakup  starts,  the  distribution  function  associated  with  this  parcel,  is  Dirac  function  at  radius 
of  the  parent  drop.  After  passage  of  time,  which  is  inversely  the  breakup  frequency,  the  new 
droplets  arise  due  to  breakup.  In  sequel,  the  droplet-radius  distribution  function  changes.  We 
suppose  that  the  new  distribution  may  be  described  according  to  solution  (9)  taken  at  vt  = 1 

with  and  (£,2^  as  parameters  of  the  model.  In  order  to  alleviate  computations,  we  can 
proceed  the  following  way.  Let  us  assume  that  once  every  breakup  time  scale,  all  outcomes  of 


breakup  in  the  given  parent  parcel  are  in  mean  (over  many  computations),  recovered  by  one 
new  parcel  that  replaces  the  parent  one.  The  radius  of  droplet  associated  with  produced  parcel 
is  sampled  from  (9).  The  new  number  of  droplets  is  computed  by  mass  conservation  from  the 
primary  parcel  to  the  secondary  one.  After  the  sampling  procedure,  the  current  time,  t , 
prescribed  for  produced  parcel  is  counted  from  zero  and  Lagrangian  tracking  is  continued  up 
to  the  moment  {vt  = 1 ) of  the  further  breakup. 

In  computations,  we  used  expressions  obtained  for  the  distribution  of  the  logarithm  of  radius. 
The  starting  distribution  for  the  logarithm  of  droplet  radius  in  j - th  primary  parcel  is 


^o/(*o)=  <?(*<>-*/) 


(10) 


Using  this  distribution  function  in  (7)  at  vt  = 1 , one  can  express  the  solution  by  the  error 
function  erf : 


r,M=r 


\ + erf\ 


m . 


(ii) 


The  product  droplet  velocity  is  computed  by  adding  to  the  primary  parcel  velocity  a velocity 
whu , which  is  randomly  distributed  in  a plane  normal  to  the  relative  velocity  vector  between 

the  parent  droplet  and  gas.  The  quantity  of  |wAu|  is  determined  by  the  mean  local  Sauter 
radius  of  parent  drops,  r, , and  the  breakup  frequency,  v : 

W = /*3,v  (12) 

IV.2  Critical  radius,  breakup  frequency 

The  critical  (or  maximum  stable)  radius  is  determined  when  disruptive  hydrodynamic  forces 
are  balanced  by  capillary  forces: 


rcr  = Wecr5l  pgur2  (13) 

where  ur  is  the  relative  between  liquid  and  gas  velocity,  8 is  the  surface  tension  coefficient, 
Wecr  is  the  critical  Weber  number,  which  can  be  taken  of  order  one  over  a large  interval  of 
Ohnesorge  numbers  l0  M.  In  the  paper,  written  in  1949  l2,  Kolmogorov  considered  the 
stretched  drop  of  insoluble  liquid  that  was  submerged  in  a turbulent  flow.  Using  the  Obuchov- 
Kolmogorov’s  scaling  law  for  the  velocity  difference  across  a size  when  the  surface  tension 
force  becomes  significant,  Kolmogorov  introduced  a critical  size  of  produced  droplets  as: 
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where  £ is  the  mean  viscous  dissipation  rate  and  p is  density  in  the  gas. 

An  estimation  of  rT  by  using  experimental  data  from 4 gives  an  enhanced  magnitude  of  rcr 
comparing  to  the  measurements.  In  order  to  account  for  the  inertia  of  liquid,  namely  for  the 
density  of  the  liquid  p , , the  expression  (14)  can  be  modified.  Estimating  the  mean  square  of 
relative  droplet-to-gas  velocity  by  mean  viscous  dissipation  and  Stokes  time  scale1314 , 


(iir2)-«-„ 


one  yields  a new  expression  for  critical  radius: 


36,/3  f WeSv 


2 { epl 


Using  the  experimental  data  from  4:  water  density,  1000  kg/m3;  gas  viscosity,  1.5xl0'5m2/s; 
gas  orifice  size,  2.1  mm;  surface  tension,  0.07  kg/s  ; gas  injection  velocity,  140  m/s  and  by 
setting  the  turbulent  gas  velocity  at  one  tenth  of  the  gas  injection  velocity,  one  gives  for 
critical  radius  3x1 0'5  m,  which  is  of  the  same  order  that  was  measured  in  4.  At  the  same  time, 
expressions  (17)  requires  a reliable  knowledge  of  viscous  dissipation  rate,  which  is  a problem 
in  the  turbulence  computation.  For  these  reasons,  the  critical  radius  is  calculated  in  this  work 
by  the  standard  expression  (13),  where  ur  = v - \p  is  calculated  by  the  mean  relative 

velocity  between  gas  and  liquid  particle,  computed  by  the  model  of  turbulent  dispersion  of 

1 3 P ( W, 

particles  . Note  that  introducing  the  turbulent  Weber  number,  Welur  = g — - - , and  using 

8 

(16),  one  may  write  for  the  critical  Weber  number: 


We  -——Re  We  \ — 

36  Pz  l L 


Assuming  that  at  scales  where  breakup  takes  place  Retur  is  of  order  of  unity  and  llur  ~ rj , one 


may  propose: 


We,  W p 


-=3.3  — 2- 
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The  choice  of  the  breakup  frequency  has  to  be  stated  from  the  physics  of  atomization.  In  this 
paper,  the  breakup  time  scale  is  taken  from  15~17: 


v = B> 


v„  - v 


A P * 


where  rr  is  the  local  Sauter  mean  radius  of  parent  drops  and  5 = 1/  V3  is  taken  from  TAB 
model  ,8. 

IV  J Choice  of  parameters  and  {^2j 


Multiplying  (8)  by  r and  integrating  over  the  entire  r - range  give  us  an  expression  for  the 
first  moment 

(r) = (rL,  expiK  fe) + °-5(£2))  t\  & 1 ) 

The  condition 

® < -{(?)  W 

provides  for  ~ < 0.  In  this  paper,  the  magnitude  for  U2)  is  supposed  to  be  proportional  to 
the  maximal  dispersion  of  radius  (£2)oc  In  1-ln— . Replacing  in  (19)  r\  by  the  local  Sauter 

r32 

mean  diameter  of  parent  drops,  one  may  assume  that 


(^2^  ~ -In—  = const  ■ In 


and  (c)  is  an  arbitrary  parameter  to  be  taken  according  to  (22)  and  (23). 


V.  EXAMPLE  OF  LAGRANGIAN  COMPUTATION  OF  THE  ATOMIZING 

SPRAY 

The  configuration  and  inlet  parameters  from  the  CORIA’s  injector19  are  used  in  the 
computation.  In  this  experiment,  the  round  jet  of  water  issues  from  the  central  tube 
( D,  = 1 .8 mm ) at  low  velocity  and  atomizes  by  a parallel  coflow  of  air  issuing  at  high  velocity 

from  an  annular  duct  ( D K = 3.4 mm ).  An  example  of  spatial  distributions  of  drops  in  the  spray, 

the  schematic  of  the  injector  and  the  evolution  in  time  of  distributions  of  droplet-size 
probability  density  function  (pdf)  at  two  cross  sections  in  the  near-nozzle  region  are  given  in 
Fig.l.  The  pdf  distributions  are  scaled  on  the  total  drops  number  crossing  the  given  section  at 
the  given  time  moment.  The  statistics  of  radius  at  (3-5)  mm  shows  mostly  the  large  unbroken 
drops  of  size  of  the  injector  orifice  that  are  accompanied  by  small  striped  droplets.  From  pdf  s 


at  (7-9)  mm,  it  is  seen  that  the  probability  to  find  drops  of  size  of  the  injector  orifice  is 
essentially  decreased  while  new  droplets  are  produced  with  radius  from  50  pm  to  250  pm. 
These  figures  and  zooming  given  in  Fig.2  show  that  a broad  spectra  of  droplet  size  is 
presented  by  computations  with  a co-existence  of  large  drops  and  small  droplets. 
Computations  are  performed  for  different  inlet  air  and  water  velocities  represented  in  19 

p U~ 

providing  for  the  different  magnitudes  of  the  parameter  J = — f-  and  the  inlet  Weber 

Piui 

number.  The  qualitative  agreement  have  been  obtained  between  measured  impact  liquid  core 
and  the  estimated  length  of  the  zone  presented  by  computed  blobs  of  size  of  the  injector 
orifice.  The  modified  numerical  code  with  the  stochastic  model  of  breakup  is  specifically 
target  on  the  computation  of  the  spray  combustion  in  the  configuration  likely  to  rocket  engine. 

VI.  CONCLUSION 

A new  sub-grid-scale  stochastic  model  of  drops  air-blast  breakup  is  presented  in  this  paper. 
The  stochastic  process  is  considered  in  the  framework  of  cascade  of  uncorrelated  breakage 
events  in  series  down  to  the  critical  size,  independently  from  the  initial  distribution  of  sizes. 
To  this  end,  the  Kolmogorov’s  discrete  model  of  particle  breakup  has  been  reproduced  in  the 
form  of  evolution  equation  for  distribution  function.  The  asymptotic  solution  of  this  equation 
has  been  applied  to  simulate  the  drop  breakup  alongside  with  Lagrangian  model  of  spray 
dynamics.  Performed  computations  of  air-blast  atomization  are  related  to  a spray  close  to  the 
rocket  engine  configuration.  A broad  spectra  of  droplet  size  is  simulated  at  each  spray 
location  with  a co-existence  of  large  drops  and  small  droplets.  The  evolution  of  shape  of  the 
droplet-size  pdf  s is  shown  at  different  sections  in  downstream  direction. 
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Fig.l  Computed  one-side  spatial  distributions  of  drops  and  PDF’s  of  size 
at  different  sections.  CORIA  GDR  injector. 
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